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We present large-scale molecular dynamics simulations of a nematic-isotropic interface in a system 
of repulsive ellipsoidal molecules, focusing in particular on the capillary wave fluctuations of the 
interfacial position. The interface anchors the nematic phase in a planar way, i.e., the director 
aligns parallel to the interface. Capillary waves in the direction parallel and perpendicular to the 
director are considered separately. We find that the spectrum is anisotropic, the amplitudes of 
capillary waves being larger in the direction perpendicular to the director. In the long wavelength 
limit, however, the spectrum becomes isotropic and compares well with the predictions of a simple 
capillary wave theory. 
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I. INTRODUCTION 

Surfaces and interfaces in liquid crystals have been the 
subject of much interest both from a fundamental point 
of view and because of their practical importance in the 
context of liquid crystal devices . 

Liquid crystals are formed by anisotropic particles. De- 
pending on temperature and density, they exhibit vari- 
ous liquid phases. Here, we shall be concerned with the 
isotropic phase (I), which is an ordinary fully symmet- 
ric fluid phase, and the nematic phase (N), where the 
fluid retains the translational symmetry in all directions 
(no positional order), but exhibits long range orienta- 
tional order. According to a very general symmetry ar- 
gument, the transition between the isotropic and the ne- 
matic phase is first order Q. Therefore there exist re- 
gions in phase space where nematic and isotropic fluids 
coexist and are separated by interfaces. 

The nematic state is characterized by a so-called ne- 
matic director, which specifies the preferred direction of 
alignment of the particles. Surfaces and interfaces break 
the isotropy of space and usually favour a certain di- 
rector orientation. This effect, called surface anchoring, 
is commonly characterized by the tilt angle 6 between 
the preferred orientation and the surface (interface) nor- 
mal, and the anchoring strength or anchoring energy. In 
particular, the cases 6 = 0° (perpendicular alignment) 
and 9 = 90° (parallel alignment) are often referred to 
as homeotropic anchoring and planar anchoring, respec- 
tively. 

Theoretically, particular attention has been given to 
Nl-interfaces in systems of hard particles. Studies of NI- 
interfaces in hard rod systems based on Onsager's ap- 
proach H have indicated that the surface free energy 
has a minimum when the rods lie parallel to the inter- 
face (|,U, suggesting that the anchoring in this system 
is planar. Similar results have been obtained in systems 
of hard ellipsoids |6) . According to the Onsager theory, 
the interfaces are flat, i.e., fluctuations of the interface 



position (capillary waves) are not present. 

In general, however, one would expect such fluctua- 
tions in fluid-fluid interfaces. The usual simple argument 
reads as follows: On very large length scales, the free en- 
ergy of a system containing an interface should be gov- 
erned by the interfacial tension 7. Let us neglect bub- 
bles and overhangs and parametrize the local position of 
a fluctuating planar interface by a single- valued function 
h(x, y). Fluctuations enlarge the interfacial area and thus 
cost the free energy (7]|| 




(assuming small distortions with \dh/dx\, \dh/dy\ <C 1). 
Since this free energy functional is quadratic, the par- 
tition function can be evaluated exactly. ^F{h} is eas- 
ily diagonalized by means of a two dimensional Fourier 
transform. One finds that the squared amplitude of fluc- 
tuations with wavevector q is on average given by 

<IMq)| 2 H^?, (2) 
79 

i.e., diverges in the long-wavelength limit q — > 0. The ar- 
gument thus not only predicts on fairly general grounds 
the existence of capillary waves at all temperatures, but 
also that they are actually quite large. 

On the other hand, it has been argued that the On- 
sager approach should be exact in the limit of infinitely 
elongated particles || . The Onsager theory assumes that 
the structure of the fluid is entirely determined by the 
second virial coefficient, i.e., by the statistical properties 
of clusters of two particles. In systems of very long, very 
thin particles, the probability that one particle interacts 
with more than one other particle approaches zero, and 
the Onsager assumption becomes exact. One concludes 
that capillary waves should vanish in this limit, in con- 
tradiction to the argument presented just above. 
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How can this be? Two explanations are conceivable. 
First, capillary waves vanish if the interfacial tension di- 
verges. However, this is unlikely here, since the coex- 
istence density of the nematic and the isotropic phase 
moves towards zero as the particle elongation is in- 
creased. Second, capillary waves might be suppressed 
by long-range interactions, which are disregarded in the 
functional (|l|) . Indeed, elastic long-range interactions are 
present in the nematic phase. We recall that the inter- 
face orients the director of the nematic liquid in a cer- 
tain direction. Long-range interface fluctuations there- 
fore lead to long-range elastic distortions of the director 
field, which are penalized. 

In order to shed more light on these issues, we have 
carried out extensive Molecular Dynamics simulations of 
an N-I interface in a system of over 100000 particles of 
elongation 15, and analyzed the capillary wave spectrum 
in detail. 

To our knowledge, there exist only few numerical stud- 
ies of N-I interfaces |Fj^MO . In one of them pJ, a 
system of molecules with length-to-width ratio k = 3 
interacting via a Gay-Berne potential jl2| was studied 
in a simulation box with a temperature inhomogeneity, 
which served to make and maintain the NI interface. The 
molecules in the nematic phase were found to align par- 
allel to the interface (planar anchoring). However, due 
to the use of the temperature inhomogeneity, the system 
is not at thermodynamic equilibrium. Capillary waves 
of the interface position fluctuations are obviously sup- 
pressed. In the other studies ||[ll], repulsive ellipsoids 
of revolution of axis ratio k = 15 were used in a sys- 
tem confined between two hard walls. The interactions 
of the walls with molecules were systematically varied 
to study the interplay between surface and interface an- 
choring. Planar alignment at the N-I interface was again 
observed. The system sizes considered were too small to 
allow for an analysis of capillary waves. 

The present work builds on Reference [jnj. We study 
Nl-interfaces in a very large system of repulsive ellipsoids 
by means of extensive computer simulations, focusing in 
particular on the capillary wave fluctuations of the in- 
terfaces. The paper is organized as follows: In the next 
section, we describe the model and some simulation de- 



tails . Section III explains our way to determine the local 
position of the interface. The results are described and 



summarized in section IV. 



T/(u l ,u J ,r u ) = |/ otherwise ^ 



with 



X = 



r - cr(Uj,Uj,fy) 



(4) 



Here, and Uj are the orientations of ellipsoids i and j, 
and Yij is the center-center vector between the ellipsoids, 
of magnitude r and direction f ^ . The distance function 
a [^3) approximates the contact distance between two 
ellipsoids and is given by 



r V 

cr(Ui,Uj-,f) = o- s jl - - 



r) 2 



(•5) 



1 - %u 4 ■ u 3 

where the parameter \ is related to the elongation k 
ai/a s through 
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Since the attractive tail of the potential eqn. (||) has been 
cut off, its range is much shorter than that of the Gay- 
Berne potential. This enables us to simulate systems 
with a large number of particles, which is essential to 
study the capillary wave effects of interest to us. For 
convenience, a s = 1 defines a unit of length, e = 1 a unit 
of energy, and we take particle mass m = 1; the particle 
moment of inertia is set at / = 50mcr 2 . 

As a preliminary run, we performed a molecular dy- 
namics (MD) simulation of a system with 7200 molecules 
at constant volume in a box geometry (L x : L y : L z ) — 
(1:1:8) with periodic boundary conditions in all direc- 
tions. The density was chosen in the coexistence region, 
p = 0.017/erf, and the temperature k^Tjt = 1. A ne- 
matic slab bounded by two interfaces parallel to the xy 
plane assembled as a result. More than 5.5 x 10 6 MD 
steps were required to equilibrate this system, where each 
step covers 0.002 time units. The last configuration was 
then reproduced 4 times in the x and y direction, which 
generated a system with N = 115200 molecules and box 
size L x = L y = 150.1958cr s = L and L z = 300.3916cr s . 
This was used as the starting point for a MD simulation 
of 4.2 x 10 6 MD steps. 



II. MODEL AND SIMULATION DETAILS 

Our system consists of idealized ellipsoidal particles 
with elongation n = ai/a s = 15 where 07 and a s are the 
length and width of the particles, respectively. The large 
value of the elongation k = 15 ensures that the order 
of the transition is fairly strong, which makes it easier 
to generate and maintain the interface. The interaction 
between two ellipsoids i and j is given by 



III. BLOCK ANALYSIS AND THE DIVIDING 
SURFACE 

In order to study the interfacial position fluctuations 
as a function of the system size, one can either perform 
simulations of various system sizes or simulate a single 
(large) system and analyze subsystems of it. In this 
study, we took the latter approach. We split our sys- 
tem of size L x L x L z into columns of block size B x B 
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FIG. 1. Illustration of our method of estimating the divid- 
ing surface. Solid line shows an example of an order param- 
eter profile for one of 16 columns with B = L/4. Dotted, 
dashed and long dashed lines show averages of this profile 
over iV avo = 15, 30, and 45 bins, respectively. We defined the 
position of the interface z- m t to be the intersection of the two 
profiles averaged with -/V avc = 15 and N me = 30. 



where Ui a is the a component of the unit vector which 
points along the axis of the i-th. molecule (a = x,y,z), 
and n is the number of molecules in the bin. The nematic 
order parameter S(z) in the bin centered at z is given by 
the largest eigenvalue of the ordering tensor. The further 
procedure is illustrated in Fig. [l]. Once the order param- 
eter profile S(z) is obtained (solid line), we compute at 
least two averaged profiles over iV avc bins, where -/V avc is 
chosen such that an averaged profile still reflects the ex- 
istence of the interface, but short range fluctuations are 
averaged out (dotted and dashed lines in Fig. lj). Finally, 
the position of the interface is estimated as the intersec- 
tion of two averaged profiles. 

The method is motivated by the following considera- 
tion: If the order parameter S(z) did not fluctuate at all 
in the bulk, the intersection would locate a "dividing sur- 
face" , where the negative order parameter excess on the 
nematic side just balances the positive order parameter 
excess on the isotropic side. It should then be indepen- 
dent of iVave, as long as L N avc /Nb ins is larger than the 
interfacial width and smaller than the total width of the 
nematic slab. This is not exactly true in an actual con- 
figuration due to the bulk fluctuations of S(z), but the 
methods still works well even for small block sizes (see 
Fig-!)- 
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FIG. 2. Same as Fig. 1 but for B = L/10. The fluctuations 
in the order parameter profile are much larger. 



and height L z . The columns are further divided into 
iVbins = 200 bins in the z direction. The local position of 
the interface Z[ nt (x,y) flij in each column is estimated 
as follows. 

First, we compute the local ordering tensor S in each 
bin (of size B x B x (L z /Ni,- ms )), defined as 



IV. RESULTS 

The fluctuations of the total director (the direction of 
the eigenvalue corresponding to the largest eigenvalue of 
the total ordering tensor) turned out to be so slow that 
the director hardly changed throughout the run. It al- 
ways pointed in the y direction. This made it convenient 
to resolve wave-vector components along and perpendic- 
ular to the director, without the need to apply any kind 
of director constraint. 
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S a /3 = — ^ 2 (3Ui a Ui/3 — Sap), 
i—1 



(7) 



We used the procedure sketched in section [II to de- 
termine the local deviations h{x,y) = Zi nt {x,y) — z int of 
the local interface position from its average for various 
block sizes B. The landscape obtained for the block size 
B = L/8 was further analyzed and Fourier transformed. 
Fourier modes h(q) are labeled by n — (n x ,n y ), where 
n x and n y are positive integers with q = -?n. 

First, we inspect the relaxation times and the correla- 
tion times of the Fourier modes. The time evolution of 
the modes n = (0, 1), (1, 0), (0, 4), and (4, 0), is shown as 
a function of time (in MD steps) in Fig. 0. From there 
we estimate the number of MD steps needed to equili- 
brate the system. For the slowest mode n = (0,1), in 
the direction parallel to the director, the equilibration 
process seems to require roughly 1.0 x 10 6 MD steps. 
Based on this information, we have discarded the initial 
1.2 x 10 6 MD steps, and collected results over the follow- 
ing 2.96 x 10 6 MD steps only. 

One also has to ensure that the total length of our sim- 
ulation run significantly exceeds the characteristic time 
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„ FIG. 5. Distribution of local interface positions for block 

FIG. 3 Fourier modes |fc(q)| for n = (0, 1), (1, 0), (0, 4), sizeg g = L/ ^ L/? and L/g The soM Unes arg Qaus _ 

and (4, 0) vs. time in units of MD steps. g j an g^_ g , 
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FIG. 4. Autocorrelation function Cuh (log scale) vs. MD 
steps for the two slowest modes n = (1,0) (normal to the 
director) and n = (0, 1) (parallel to the director). Dashed 
lines are linear fits to give the correlation time. 



scale of the slowest capillary wave mode. In order to 
check this, we have computed the autocorrelation func- 
tions of the squared amplitude of the Fourier modes /i(q) 



C hh (t) = 



<|ft(t)| 2 |fc(0)| a >- 



|2\2 



-t/T 



(8) 



(N 4 }-(N 2 ) 2 

They are shown in Fig. |I] for the slowest Fourier modes 
in the direction normal to the director n = (1,0) and 
parallel to the director n = (0,1). The correlation time is 
estimated from the slope of the fitted line as r ss 1.0 x 10 5 
MD steps for the (1, 0)-mode and r » 3.0 x 10 5 MD steps 
for the (0, l)-mode. This is only by a factor of 10 smaller 



than the length of the total simulation run. Hence the 
statistical error of our results is quite large. 

We can now proceed to compare our results with the 
predictions of the capillary wave theory. 

First, we consider the distribution of h(x, y). From the 
free energy functional ([!]), one derives the exact predic- 
tion 111 
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where 
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~ 4^2 



dq(\h(q)\ 2 
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The lower and upper cutoffs q ma ,x — 2tt/B and 
<7min = 2tt/L come into play, because the integral 
J dq{\h(q)\ 2 ) ~ J dq/q diverges as q — > and q — > oo. 
The actual local height distribution P(h) obtained from 
our simulations is plotted for several block sizes in Fig. |^. 
It can be fitted nicely by a Gaussian distribution @. 

Next, we study the width u of the average order pa- 
rameter profile as a function of the block size. According 
to the capillary wave theory, the capillary wave fluctua- 
tions broaden it for large blocks B according to 
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(11) 



Here the cutoff wavevectors are given by q mm — 2ir/B 
and g max = 27r/ao, where ao is a microscopic length which 
need not be specified here. 
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FIG. 6. Order parameter profiles for block sizes B = L, 
L/2, L/4, I//6 and L/8. Sn and Si are the average values 
of the order parameter in the nematic and isotropic phase, 
respectively. 
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FIG. 7. Squared interfacial width uj 2 vs. block size B. Solid 
line is a fit to eq.(ll). 



The broadening effect is demonstrated for different 
block sizes B in Fig. ^. The interfacial width can be es- 
timated by fitting order parameter profiles such as those 
shown in Fig. |^ to £cm/i-profiles, 



S(z) = - (S n + Si) + - (S N - Si) tanh 



Zint 



(12) 



where the parameters Sn and Si are the values of the or- 
der parameter in the bulk nematic and isotopic phases. 
(Note that Si is nonzero due to the finite number of par- 
ticles in a bin). Fig. ^ plots the squared interfacial width 
u) 2 as a function of block size B. One clearly observes the 



logarithmic increase of lu 2 with block size B predicted by 
eqn. (|j"l|). From the slope of the line one can estimate 
the interfacial tension, 7 = 0.016 ± 0.002k B T/a 2 . 

For comparison, we have also determined the interfa- 
cial tension from the anisotropy of the pressure tensor, 
making use of the relation fl(|llj 
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dz [P N (z) - P T (z)], 



(13) 



where Pn and Pt are the normal and transverse pressure 
tensor components. Here, we consider particles with pair 
interactions , eqn (||) and systems with two interfaces 
in the (xy)-plane. Neglecting finite size effects and inter- 
actions between the interfaces, eqn. ( |l3|) thus reads 
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(14) 



The simulation data yield 7 = 0.0093 ± 0.003a 2 /k B T. 
This value agrees with earlier results obtained in smaller 
systems 17|. It is of the same order, but smaller than 
the value estimated from the interfacial broadening. The 
quantitative difference possibly stems from attractive in- 
teractions between the two interfaces of the nematic slab. 
More simulations, in which the thickness of the slab is 
varied systematically, would be needed to elucidate this 
point. 

Finally, we turn to the analysis of the capillary wave 
spectrum (|/i(q)| 2 ), which is predicted to be inversely pro- 
portional to q 2 according to eqn (j^) . We study separately 
the q-direction parallel (q oc (0,1)) and perpendicular 
(q oc (1,0)) to the director. The inverse of (|/i(q)| 2 ) 
in these two directions is shown in Fig. | as a func- 
tion of the squared wave vectors. In the long wavelength 
limit q — > 0, the spectrum appears to be isotropic, and 
l/(\h(q)\ 2 } approaches zero in agreement with the capil- 
lary wave theory. The initial slope is consistent with eqn. 
(||) if one uses 7 = 0.016fcBP/c 2 , and still roughly com- 
patible if one uses 7 = 0.0093kBT/a 2 . As q increases, 
the capillary wave spectrum becomes anisotropic. As 
one might expect intuitively, the amplitudes (|/i(q)| 2 ) are 
larger for capillary waves in the direction perpendicular 
to the director. 

Deviations of capillary wave spectra from the straight 
line predicted by the simple capillary wave theory (^) are 
often discussed in terms of higher order terms in the cap- 
illary wave Hamiltonian (Q), i.e., terms proportional to 
squares of higher order derivatives of h{x,y). For exam- 
ple, including the next-to-leading term leads to a predic- 
tion of the form l/(|ft,(q)| 2 ) oc 7^ 2 + Sq 4 , where 7 is the 
interfacial tension and 6 is a bending rigidity. One might 
intuitively expect that the Nl-interface has bending stiff- 
ness, based on the argument that the elastic interactions 
should penalize interfacial bending. 

However, Fig. || indicates that the sign of the "bending 
energy" is negative. The naive argument sketched above 
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FIG. 8. Inverse of the mean-squared Fourier components of 
the interface position l/(|/i(q)| 2 ) vs square of the wave- vector, 
n 2 = (qL/2ir) 2 for q parallel (circle) and normal (square) to 
the director. Dashed curves are guides to eye. Solid lines 
correspond to the prediction of the simple capillary wave 
theory (eqn (2)) with 7 = 0.016cr 2 /fcsT (upper line), and 
7 = 0.0093a 2 /k B T (lower line). 



is thus clearly wrong. On the contrary, our unexpected 
result suggests that the elastic interactions influence the 
capillary waves on the largest length scale: The inter- 
face is rougher on short length scales than one would 
expect from the long wavelength fluctuations. This ob- 
servation is consistent with our speculations in the intro- 
duction, that the elastic interactions may be responsible 
for the suppression of capillary wave fluctuations of the 
Nl-interface in the Onsager-limit of infinitely elongated 
particles. 

We note that this is not the first observation of a neg- 
ative bending rigidity at a fluid-fluid interface. A similar 
phenomenon has been predicted theoretically for liquid- 
vapour interfaces [|18|J20] and verified experimentally by 
C. Fradin et al [|21|| . The negative bending rigidity is 
attributed in their case to long range van der Waals in- 
teractions fl9[| . Indications of a negative bending rigidity 
at a liquid-liquid interface have also been found in Molec- 
ular dynamics simulations of simple liquids by Stecki and 
Toxvaerd ||,|§. 

To summarize, we have studied the Nl-interface by 
large scale molecular dynamics simulation of a system 
of repulsive ellipsoidal molecules. We find that the stan- 
dard capillary wave theory explains most of our results 
very well. Discrepancies are encountered when looking 
at the amplitudes of capillary wave modes at large wave 
vectors, i.e., on small length scales. We find that they 
are smaller than expected and anisotropic. 
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